home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 138_01 / rkst1.c < prev    next >
Text File  |  1985-03-09  |  4KB  |  175 lines
  1. /*
  2.     program:    rkst1.c
  3.     created:    03-Nov-83
  4.     by:        A. Skjellum
  5.  
  6.     modified:    14-Nov-83
  7.     by:        M. J. Roberts
  8.  
  9.     Copyright 1983, 1984 (c) California Institute of Technology.
  10.     All rights Reserved.  This program may be freely distributed
  11.     for all non-commercial purposes but may not be sold.
  12.  
  13.     purpose:    illustrate the use of rk4n program
  14.  
  15.     uses:        rk4n() (rks.c)
  16.             
  17.     summary:
  18.         integrates the differential equation:
  19.  
  20.         y'(t) + y(t) = t + 1
  21.         y(0) = 5.0.
  22.  
  23.         for which the exact solution:
  24.  
  25.         y(t) = t + 5exp(-t) is known.
  26.     
  27.         Integrates the same equation as rktest1 
  28.         but using the more general equation solver, 
  29.         rk4n().  This run is exactly the same as for 
  30.         rktest1, except that, rather than trying to solve 
  31.         exactly the same equation, we will solve a 
  32.         "system" of one differential equation.
  33.  
  34. */
  35.  
  36. /* constants */
  37.  
  38. #define SYSIZE    1    /* number of functions in system */
  39. #define    YZERO    5.0    /* initial value for y */
  40. #define    TSTART    0.0    /* starting time for integration */
  41. #define    TEND    10.0    /* ending time for integration */
  42. #define    STEPS    50    /* 50 steps in integration */
  43.  
  44. /* variables external to all functions */
  45.  
  46. double wvalue[STEPS+1][SYSIZE];
  47. double yarray[STEPS+1][SYSIZE];
  48. double tarray[STEPS];
  49.     /* integrated solution stored here */
  50.  
  51.  
  52. /* subroutines: */
  53.  
  54. /* exact(): returns exact solution value, given t */
  55.  
  56. double exact(t)
  57. double t;
  58. {
  59.     extern double exp();    /* exponential function */
  60.  
  61.     return((t + YZERO*exp(-t)));
  62. }
  63.  
  64. /* fn(j,i,t,y): return f(t,y) given t,y values */
  65.  
  66. double fn(j,i,t,y)
  67. int j;
  68. int i;
  69. double t;
  70. double (*y)();
  71. {
  72.     double a,b;        /*  temporary storage space  */
  73.  
  74.     /*
  75.         differential equation is y' + y = t + 1
  76.  
  77.         therefore, y' = t + 1 - y.
  78.  
  79.     */
  80.  
  81.     a = (*y)(0,j,i);    /*  calculate function
  82.                 Note that the ZERO was passed so as to
  83.                 allow the function to know which argument
  84.                 we are talking about -- in this case,
  85.                 we only need one argument evaluated, so
  86.                 pass it 0 to indicate the first (zeroeth,
  87.                 actually) argument is to be calculated.  */
  88.  
  89.     b = t + 1.0 - a;    /*  and figure out the rest of it  */
  90.     return(b);
  91.  
  92. }
  93.  
  94.  
  95.  
  96. /* store():  the routine to store away the W values for later reference */
  97.  
  98. double store(row, col, value)
  99. int row, col;        /*  location to store the value  */
  100. double value;        /*  the actual value to store  */
  101. {
  102.     wvalue[row][col]=value;
  103.     return (value);
  104. }
  105.  
  106. /* source():  return the  W value referenced by input parameters  */
  107.  
  108. double source(row,col)
  109. int row, col;        /*  location to look up */
  110.  
  111. {
  112.     return (wvalue[row][col]);
  113. }
  114.  
  115. /* solutn(): print solution step at console */
  116.  
  117.  
  118. solutn(j,i)
  119. int j,i;    /* element numbers  */
  120. {
  121.     double time;
  122.     double ex;
  123.     double approx;
  124.  
  125.     time   = tarray[j];
  126.     ex     = exact(time);
  127.     approx = source(j,i);
  128.  
  129.     printf("t = %7.3e, y = %7.3e, y_exact = %7.3e, diff = %7.3e\n",
  130.         time,approx,ex,approx - ex);
  131. }
  132.  
  133. /* main program: */
  134.  
  135. main()
  136. {
  137.     /* external declarations */
  138.  
  139.     double store(), source();
  140.  
  141.     double fn();    /* ensure that this is typed as double */
  142.  
  143.     /* local variables: */
  144.  
  145.     register int i,j;
  146.  
  147.  
  148.     double init[1];        /* initial condition matrix */
  149.  
  150.     /* begin code: */
  151.  
  152.      printf("\n\nrktest1.c    as of 03-Nov-83\n\n");
  153.     printf("Integrates: y' + y = 1 + t     for\n\n");
  154.     printf("t = %7.3e to %7.3e, with %u steps\n\n",
  155.         TSTART,TEND,STEPS);
  156.  
  157.     /*
  158.         integrate the answer from t = 0 to t = 10 sec 
  159.         STEPS points.
  160.     */
  161.  
  162.     init[0] = YZERO;    /* set up initial condition matrix */
  163.  
  164.     rk4n(fn,source,store,SYSIZE,TSTART,TEND,STEPS,init,tarray,yarray);
  165.             /* compute the answers for 1 function */
  166.  
  167.     /*  Print out solution  */
  168.  
  169.     for(j=0;j<STEPS;j++)    /* print solution */
  170.         solutn(j,0);
  171.  
  172.  
  173.     printf("\n\nEnd of execution\n\n");
  174. }
  175.